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Abstract. The localization length L2 of two interacting particles in a one-dimensional disordered system is 
studied for very large system sizes by two efficient and accurate variants of the Green function method. The 
numerical results (at the band center) can be well described by the functional form L2 = Li [0.5 + c(U) Li] 
where L\ is the one-particle localization length and the coefficient c(U) « 0.074 |f/|/(l + \U\) depends on 
the strength U of the on-site Hubbard interaction. The Breit-Wigner width or equivalently the (inverse) life 
time of non-interacting pair states is analytically calculated for small disorder and taking into account the 
energy dependence of the one-particle localization length. This provides a consistent theoretical explanation 
of the numerically found [/-dependence of c(U). 

PACS. 72.15.Rn Quantum localization - 71.30.+h Metal-insulator transitions and other electronic tran- 
sitions - 72.10.Bg General formulation of transport theory - 71.55.Jv Disordered structures; amorphous 
and glassy solids 



1 Introduction 

The quantum eigenstates of non-interacting particles in a 
random potential are localized if the fluctuations of the po- 
tential (the disorder strength) are sufficiently strong (for 
a review see jp). This phenomenon of Anderson localiza- 
tion is particularly well understood for one-dimensional or 
quasi one-dimensional geometries where localization even 
persists for arbitrarily small disorder strength. For this 

■ case efficient numerical methods [Q and also powerful an- 
alytical theories in terms of the supersymmetric non-linear 
er-model and the Fokker-Planck approach for the 

transfer matri x B pl are available. 

Dorokhov ||and, recently, Shepelyansky || consid- 
ered the case of two interacting particles (TIP) moving in 
a one-dimensional random potential for which they pre- 
dicted a strong enhancement of the localization length for 

' the pair-states due to the interaction. While the results 
of Dorokhov are only valid for the case of a strongly at- 
tractive interaction confining both particles together, She- 
pelyansky also considered a local, attractive or repulsive 
Hubbard interaction. He claimed that among many states 
with both particles being localized far away, there are a 
few pairs-states where the typical distance between the 
particles is of the order of the one-electron localization 
length L\ and the center of mass coordinate is character- 
ized by a pair-localization length L2 S> L\. Mapping the 
original problem on a random band-matrix model super- 
imposed with large diagonal elements, he found for the 
case L\ ^S> 1 (length measured in units of the lattice spac- 



ing) 



L, 



El 
t 2 



(i) 



where U is the strength of the Hubbard interaction and t 
is the value of the hopping matrix element. The disorder 
strength enters through the value of L\ (see below). This 
estimate has been confirmed by Imry |Q using a different 
argument based on the Thouless scaling block picture || . 
A crucial role is played here by the spreading width F (also 
called Breit-Wigner width) which is the energy scale over 
which unperturbed states are mixed due to the interac- 
tion. Imry identified this energy scale with a generalized 
Thouless energy defined as the sensitivity of the energy 
levels with respect to a change of the boundary conditions 
for a finite block of size L\. The pair-localization length is 
then obtained by scaling theory as L2/L1 = F j Ai where 
A2 ~ t/L\ is the two-particle level spacing in the finite 
block. Using an ergodic hypothesis for the one-particle 
eigenfunctions, one can estimate the spreading width by 



Fermi's golden rule r ~ U / (t L\) |]6|,0,P| reproducing Eq. 

First numerical studies in terms of finite size trans- 
fer matrix calculations |To| and exact diagonalization [ pd| 
confirmed the str ong enhancement due to the interaction. 
Borgonovi et al. |12j showed that the enhancement effect 
also appears in a related model of two interacting kicked 
rotors for which it is possible to determine directly the 
quantum time-evolution. Von Oppen et al. |d| introduced 
an efficient method to calculate the two-particle Green 
function and based on their numerical results they pro- 
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posed the scaling relation L2/L1 ss 0.5 + 0.054 \U\ L\ with 
a linear dependence on U contradicting the estimate (|l|). 
This behavior was explained by Jacquod et al. |Q who 
calculated analytically the spreading width for the limit 
of vanishing disorder. Resumming an infinite series of dia- 
grams they obtained for energies close to the band center 
r/A 2 ~ Li \U\/y/t 2 +JU/4) 2 . Therefore the physical ar- 
guments of Refs. did not contradict the results of 
||l3f but the application of Fermi's golden rule correspond- 
ing to lowest order perturbation theory and the ergodic 
hypothesis appeared to be insufficient to determine r. 

It is worth mentioning that the topic also inspired 
considerable progress in the understanding ||^,|^[l7]] of 
the random band matrix model with superimposed diag- 
onal originally introduced and used by Shepelyansky B. 
Furthermore, in Ref. a sophisticated random matrix 
model was proposed which works for arbitrary space di- 
mension and takes properly both particle coordinates (rel- 
ative and center of mass coordinate) into account. This 
model can be mapped onto an effective supermatrix non- 
linear cr-modcl and it is thus possible to explain features 
like a logarithmically suppressed diffusion or a logarith- 
mically increasing pair size Q found previously by Bor- 
gonovi et al. |l^] for the pair-diffusive regime in d > 2 
dimensions where all states are delocalized. Subsequent 
work was concerned with the role of the level statistics 
JRH and, very recently, with the fractal structure of the 
coupling matrix elements due to the interaction J2(J . 

Despite the available evidence in favor for the enhance- 
ment effect the general situation is still not really clear, 
due to different proposals for the dependence of L 2 on W 
and U |2(],|2l],^2| and the claim of Romcr et al. that the 
effect completely vanishes in the limit of infinite system 
This claim, which was contested in 



size 



23 



is based 

on the finite size extrapolation of the localization length 
calculated by a transfer matrix method for finite square 
samples being put together to an infinite strip. 

In this work, we present numerical results (section |J) 
based on an exact and efficient variant of the Green func- 
tion method introduced in Ref. jl3| that allows to treat 
rather large system sizes, i. e. N > 1000. This is indeed 
important for small disorder values in order to perform an 
accurate finite size extrapolation. We furthermore present 
a second variant which consists of the recursive Green 
function technique applied to an effective band matrix 
Hamiltonian as considered in jl3|]. In this approach one 
can indeed take the limit A — > 00 and the results we find 
are consistent with those of the finite size extrapolation of 
the first variant. The issue of an accurate variant of the 
Green function method is actually of considerable interest 



since Romer et al. [ 
of von Oppen et a 



I had questioned the original results 
due to a certain approximation ap- 
plied in the original approach of Ref. (r| . We find in our 
calculations qualitative agreement with those results con- 
cerning the strong enhancement of the localization length 
Li and the dependence on L\. However, we find neverthe- 
less a quantitative difference concerning the dependence 
on the interaction strength U which is only linear for suf- 
ficiently small U . To understand this, we reconsider the 



issue of the determination of the Breit-Wigner width r 
for small disorder (section |^) . Improving the r estimate 
of Ref. 1 14 1, we can indeed explain the modified U depen- 
dence. 

Very recently, we learned of related relevant work [ p5| 
pBf in which the TIP Green function was evaluated by a 
decimation method for system sizes up to A = 251 p5| or 
N = 300 Eft 



2 Numerical Green function approach 

We consider two particles in a disordered system interact- 
ing via a local Hubbard-interaction and characterized by 
the following tightbinding Hamiltonian, 



H = -*y"l( \x + l,y>< x,y\ + \x,y + 1 >< x,y\ + h.c 
+ J2(£(x)+e(y) + U5 x , y )\x,y><x,y\ . (2) 



x and y denote the positions of the first and the second 
particle, respectively, t is the strength of nearest neighbor 
coupling matrix element which we put to unity in the fol- 
lowing and U is the value of the on-site interaction. The 
disorder energies are random, i. e. e(x) S [-W/2, W/2] 
with W being the disorder strength. At vanishing interac- 
tion U = the one-particle eigenstates (at a one-particle 
energy E = 0) are localized with the localization length: 
Li re 105/ W 2 0. In this work we do not discuss the par- 
ticular effects of symmetric or anti-symmetric two-particle 
states (Bosons or Fermions). The on-site Hubbard inter- 
action only acts on the subspace of symmetric states and 
our results apply therefore to the case of Bosons. However, 
for the actual calculations and the representation, we find 
it more convenient to keep all states. 

To determine the two-particle localization length, we 
consider as in Ref. jl3| the two-particle Green function. 
Since we are interested in the coherent propagation the 
particles being close, we determine only the Green func- 
tion matrix elements of doubly occupied states \xx >, 



g xy =< xx\ (E — H) 1 \yy > 



(3) 



A priori, for a finite system of size A, the matrix inverse 
in (||) has to be evaluated for a A 2 x A^ 2 matrix. Fortu- 
nately, von Oppen et al. have shown that this prob- 
lem can be reduced to an effective Green function on an 
A-dimensional space because the interaction operator is 
proportional to the projector on the space of doubly oc- 
cupied states. The matrix g in (j^) can be calculated [|l3] 
from an A x A-matrix equation 



1 



3o 



1-90 U 



where g = g 



(7=0 



(4) 



The matrix go is given in terms of the one-particle eigen- 
states tp a {x) and the one-particle energies E a via, 

(go)xy = y2,y a {x)yi3{x) - — — — ^r^0{y)<Poc{y) ■ 

(5) 
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The two-particle localization length L 2 
the exponential decay g xo ,x ~ exp(— |x 
sponding to 



is determined by 
— xq\/L 2 ) corre- 



-3- = - lira i ( In 

L 2 JV—oo AT 



9x,x+N 



9x,x 



(6) 



The ensemble average is performed over different disorder 
realizations and for practical purposes also over some ini- 
tial sites x close to one boundary. The extra denominator 
9x, x hi @ is not relevant in the limit TV — * 00 but pro- 
vides a considerable improvement if (^) is evaluated for 
finite N. For vanishing interaction, we expect according 
to (|) L 2 {U = 0) » Zo/2 |K§. 

In Rcf. |L3||, Eq. (Q) was evaluated for a finite system 
by employing two approximations. First, von Oppen et al. 
omitted the first factor go and, second, they did not evalu- 
ate the full matrix go but only a sufficiently large band on 
which they applied the recursive Green function technique 
p7|,p8| for the matrix inverse in (|J). Since go is indeed a 
band matrix of width ~ L\ both approximations seem to 
be well justified provided L\ < L 2 /2. However, the valid- 
ity of the corresponding results was seriously questioned 
by Romer et al. |23| due to these approximations and, fur- 
thermore, the limit of very small (vanishing) interaction 
cannot accurately be studied within this approach. 

We have evaluated (^) exactly without any approxima- 
tions. For this we note two important points concerning 
the numerical precision and the efficiency. First, the mul- 
tiplication of the band matrix go with the matrix inverse 
in (|4|) requires that the relative error of the exponentially 
small matrix elements of go far away from the diagonal 
is small [^9). Otherwise, Eq. (0) provides incorrect re- 
sults for L 2 . This in turn requires that the exponential 
tails of the ip a (x) are accurate over the whole length scale 
x = 1, . . . , N. According to this, we have determined the 
one-electron eigenstates by the method of inverse vector 
iteration |30) which provides the required accuracy by suf- 
ficiently increasing the number of iterations. The second 
point concerns the efficiency. Here the matrix inverse in 
(||), which costs of the order of N 3 operations, is actually 
not the limiting factor. This is due to the necessary eval- 
uation of go- The naive application of (^|) already costs 
of the order of N 4 operations. Even though this number 
can be reduced by exploiting the exponential decay of the 
(f a (x) this does not yield any significant improvement for 
small disorder values when L\ ~ 50 — 100. Fortunately, 
it is possible to determine go exactly with N 3 operations. 
For this we rewrite ([5]) in the form 



(go)xy = Va{x) G£j(E - E a ) <pp{y) 



(7) 
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Fig. 1. Finite size two-particle localization length L2(N) for 
U = 1.5, E = and W = 1.0 versus 1/N (100 < N < 1400). 
The full curve corresponds to the fit L^ 1 (N) = L^ 1 + C/N 
and the dash line corresponds to the fit L2 (N) — L2 + C /N. 
The insert shows 1/L2(N) versus 1/N with the corresponding 
linear fit. 



determine the limit L 2 = lirri/v^oo L 2 (N). For this, we 
have calculated the ensemble averaged inverse localization 
length L^iN) for different system sizes using the exact 
projected Green function g as given in (^). The limit for 
N — > 00 has been determined by the linear fit in 1/N of 
the inverse localization length: 



1 



L 2 {N) 



1 

T 2 



c_ 

N 



(8) 



where Gxy (E) is the one-particle Green function at energy 
E that can efficiently be determined by iV 2 operations due 
to the tridiagonal form of the one-electron Hamiltonian. 
Since this has to be done for N different energies E — E a , 
Eq. (Q) provides an algorithm with only ./V 3 operations. 

We have used two variants of the Greens function method 
The first is based on a finite size extrapolation (FSE) to 



This ansatz for the finite size extrapolation is highly sug- 
gestive from Eq. (^). Assuming, that the typical value of 
9x,x+n does not depend on N if N -C L 2 , we see that (||) 
reproduces both limits N -C L 2 and N ^> L 2 . 

The quality of the fit is indeed confirmed by the ex- 
plicit numerical values as can be seen in Fig. [l| For W — 
1.0 (Li = 105), E = 0, U = 1.5, and 100 < N < 1400, 
the fit (||) works quite well while the direct fit L 2 (N) = 
L 2 + C/N is very poor for the considered range of N val- 
ues. Of course, for L 2 <C N both extrapolation schemes 
are equivalent. However, for the case of Fig. @, L 2 and N 
are comparable and the choice of the correct method is 
crucial. 

To obtain an independent verification of the extrap- 
olation scheme (||), we have also used a second method 
which permits to determine L 2 directly for quasi infinite 
systems. For this, as in Ref. we have omitted the 

first factor in (Q) and replaced in the denominator of the 
second factor the matrix go by a block-tridiagonal ma- 
trix of block size b. This matrix can be viewed as an 
effective Hamiltonian |13|| and it is possible to perform 



the matrix inverse of 
(RGF) technique & 



) by the recursive Green function 
In contrast to von Oppen et al., 
we have applied this method for quasi infinite system size 
as described in Ref. p7| . The limit (^) can be directly 
evaluated due to the self averaging behavior of the in- 
verse localization length. The non-vanishing blocks of go 
are determined in a local approximation, i.e. to calculate 
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Fig. 2. Localization length Z/2(&) (on a logarithmic scale) for 
U — 1 and £7 = obtained by the second method [recursive 
Green function method applied on (^)] versus the inverse block 
size 1/b (40 < b < 200). The dashed lines correspond to the 
values of Li obtained by the first method (finite size extrapo- 
lation) . 



(3o)xy for x, y G {x - b + 1, . . . ,x + b} we apply Eq. 
(ffl) for a finite system of size N c > 2b containing the sites 
x £ {x - N c /2 + 1, . . . , x + N c /2}. This works if N c and 
b are sufficiently large compared to L\ because the matrix 
elements (go)xy do not depend on N c if the sites x,y are 
far away from the boundaries. For each iteration step of 
the recursive Green function procedure the disorder re- 
alization on the N c sites is shifted by b sites and go has 
to be recalculated by ([?]). We have chosen N c — 3b and 
calculated L,2(b) for different values of the block sizes in 
the interval 40 < b < 200. Here, we expect a much faster 
(exponential) convergence of L 2 as the ratio bjL\ becomes 
larger than unity. In Fig. [|, we compare for U — 1 and 
different disorder values, 1.25 < W < 3.00, the localiza- 
tion length L 2 {b) with the values obtained from the FSE 
method. 

We find overall agreement between both methods and 
L 2 (b) indeed coincides with L^ for b > b c 5 L\. For small 
disorder values it is quite difficult to arrive at this regime. 
For b < b c the values of £2(6) are typically larger than 
the values obtained by the first method. To our knowl- 
edge, the approach described above is indeed the first 
method to determine directly the TIP localization length 
for quasi infinite system size without the side effects of a 
bag-interaction This is possible, because the cutoff 

is applied on an effective Hamiltonian and the neglected 
matrix elements are indeed exponentially small if the block 
size b is sufficiently large. The results shown in Fig. || pro- 
vide therefore an additional confirmation of the validity 
of the above discussed extrapolation scheme. 

For a systematic study of the dependence on W and 
U, we used the first variant based on the FSE scheme (||) 
which appears to be more efficient, especially for large val- 
ues of L\. For the scope of this paper, we studied the band 
center E = where the localization properties are sym- 
metric with respect to the sign of U. We considered for the 
disorder values 1.0 < W < 7.0 and interaction strengths 
0.0 < U < 2.0 at least system sizes up to N = 500 and 
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Fig. 3. Enhancement factor L2/L1 as a function of L\ — 
105/W 2 . L2 is extrapolated as in ^ using system sizes up 
to iV < 1000. The considered disorder values are W = 
1, 1.25, 1.5, 1.75, 2, 2.5, 3, 4, 5, 6, 7. 



for W < 1.75 even sizes up to N = 1000. (For W = 1.0 
and U — 1.5, 2.0, we have also calculated two data-points 
with N = 1400.) Most of the data points (for the finite 
size values Li2{N)) were determined with a relative error 
smaller than 2 %. For the largest system sizes and smallest 
disorder values the relative error is 3 — 3.5 %. To verify the 
scaling relation L9/L1 w 0.5 + 0.054 suggested by 

von Oppen et al. |13|, we show in Fig. || the ratio L%/L\ 
as a function of L\ where L2 has been obtained by FSE 
from L 2 (N). 

The linear behavior in L\ is qualitatively indeed con- 
firmed but the errors for the smaller disorder values do not 
allow to exclude a behavior of the type (L2/L1 — 0.5) oc L" 
with a < 1. A corresponding fit indeed gives a ~ 0.9 but 
this depends also on the chosen offset 0.5 (the fit with the 
offset 0.55 gives a w 1.0). We mention that the slight devi- 
ations from the linear behavior can also be well described 
by an ansatz of the type (L2/L1 — 0.5) oc Li/ln(C Li) 
suggested by Borgonovi et al. p2[ . However, since the 
precision of the data does not permit to distinguish sig- 
nificantly between this and the linear behavior, we do 
not enter into more details here. For [7 = 0, we confirm 
the previous observation |^,^5| of a slight enhancement 
L 2 (U = 0)/Li « 0.5 — 0.7 which is presumably due to the 
energy average in Eq. (||) |pi^ , [25| . 

In Fig. [|, we also show the dependence of L2 on the 
disorder strength W (for U = 1 and E = 0). Previously, 
Song et al. j22| found a behavior L2 oc W~ 2 9 and as we 
can see the overall slope in Fig. [| is indeed comparable 
to this behavior. However, we find for small and large W 
values significant deviations due to the curvature in the 
curve of ln(L2) versus ln(M^). This is due to the constant 
term in the above scaling relation. Actually, the data can 
be extremely well fitted by L 2 « Li (0.55 + 0.038 Li) for 
the whole interval 1.0 < W < 7.0 and one indeed finds the 
asymptotic behavior L 2 oc W~ A for small W (L\ 3> 10) 
and L 2 oc W~ 2 for larger values of W (1 < L\ < 10). 

To extract the dependence on U, we determined the 
slope c(U) in the linear fit L2/L1 = a + c(U) L\ which 
is compared in Fig. || with the numerical data. The slope 
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Fig. 4. Dependence of L2 on W for U = 1 (double log-scale). 
The data-points are the same as in Fig. ^j. The full curve is 
the fit L 2 = Li (0.55 + 0.038 Li) with Li = 105/ W 2 . The two 
dashed lines correspond to the limiting cases L 2 = 0.038 L? oc 
W~ 4 (or L 2 = 0.55 Li oc W~ 2 ) for Li > 10 (Li < 10). The 
dotted line corresponds the behavior oc W~ 3 for comparison. 



c(U) itself is shown in the insert of Fig. || as a function 
of U. Apparently, the U dependence is not linear for the 
whole interval 0.0 < U < 2.0. This linear behavior was 
observed by v. Oppen et al. for U < 1.0 where the 
discrepancy is still quite moderate. At U = 1.0 their values 
are about 40% larger than ours. We believe that this is due 
to finite size effects and the applied approximations. 

Also the estimate c(U) oc \U\/i/l + {U/A) 2 based on 
the analytical calculation of the Breit-Wigner width for 
W = E3 only agrees with the numerical data for suffi- 
ciently small U . In the next section, we will try to explain 
this disagreement and reconsider the determination of the 
Breit-Wigner width. 



3 Breit-Wigner Width and U dependence 

The derealization effect of two interacting particles is re- 
lated to the finite life time r = T -1 of the product states 
I a/3 > with < xix 2 \af3 >= <p a (xi) <f>p(x?) The 
interaction gives rise to transitions \a(3 >— ► \y6 > which 
can be viewed as random hops of typical size L%. For short 
time scales when quantum interference effects can be ne- 
glected one therefore obtains a diffusive dynamics |^,^| 
with the diffusion constant D ~ L^t^ 1 = L\T '. Following 
a general argument developed in Ref. J3l| and applied to 
the TIP case in Refs. |6|,[|, one can estimate the localiza- 
tion length due to quantum interference effects. According 
to this the classical diffusive behavior is only valid for time 
scales smaller than the Heisenberg time corresponding to 
a wave packet of width y/Dt, i. e. 



(9) 



t < t H {t) = is cS VDt 



Here i/ c g is the density of states per length such that 
(L^eff) -1 is the level spacing in a block of size L. For the 
TIP, we have v c g f=a v 2 L\ with v 2 (E) being the energy 



Fig. 5. L2/L1 versus the scaling parameter c(U) L\ for the 
data of Fig. || (points) where the quantity c(U) is given by Eq. 
Q with A = 0.028 and k c = 0.25 (see text). The symbols for 
the different U values are the same as in Fig. |3[ The dashed line 
corresponds to L2/L1 = 0.5 + c(U) L\. The insert shows c(U) 
as a function of U. The data points are the values obtained as 
the slope in the fit L2/L1 = a + c(U) L\. The full line is the 
analytical expression ( J23| ) (same values for A and k c as above), 
the dotted line is the linear dependence 0.054 U found in Ref. 

13| and the dashed line is the dependence suggested in Ref. 

in- 



dependent two-particle density of states. ^(E) is propor- 
tional to the inverse bandwidth 1/(8 + 2W) ~ 1 with a 
logarithmic singularity at the band center (for W = 0). 
The density u e g corresponds to the number of well cou- 
pled pair states with the same center of mass coordinate 
but with different relative coordinates (up to a maximal 
value ~ L\). 

At t sa t c = VcffD, when the condition (g) ceases to 
be valid, the discrete energy spectrum can be resolved 
and localization with a localization length L2 ~ \/Dt c ~ 
v e s D sets in. This general relation for quasi Id systems 
has also been obtained in a more rigorous way using the 
supersymmetric non-linear <7-modcl [^|. In view of this, 
the Breit-Wigner width J 1 and Li are related by 



L 2 — v e s D ~ L 1 ru 2 . 



(10) 



Jacquod et al. determined T using diagrammatic per- 
turbation theory (in U) for the case of vanishing disorder, 
W = 0, and finite system size N. They argued that one 
obtains a good estimate of T for finite W by replacing 
TV with L\ due to the ballistic dynamics in Id for length 
scales up to L%. For energies close to the band center and 
moderate interaction strengths U < 1, their result reads 



\U\ 



(11) 



As we can see in the insert of Fig. 5 our numerical data 
agrees with this behavior only for very small values of U. 
We therefore feel that it is justified to reconsider the issue 
of the Breit-Wigner width which can be calculated from 
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the energy dependent local density of states 

p a a(E) = --Im < a/3\ (E + iO- H)- 1 \a(3 > , (12) 

with | a(3 > as above. Using Schur's formula to perform 
the matrix inverse, we rewrite (O) as 



(13) 



p a p(E) = --Im E + iO — E a — Ep - 



< af3\ U \a0 > +r a p/2\ 

0) + ir {1) 

P < 11 af3 

= -2 < a/3\ U (E + iO - H) U \a(3 > 



r a - r (0) + ir {1) 

1 <*0 — 1 a/3 + 1 1 a/3 



(14) 



Here U — U \ xx xx \ ^ s the interaction operator 
and H = PHP with the projector P = 1 — \a>{3 >< a(3\. 
Eq. ( [l3] ) corresponds to the Lorentzian or Breit-Wigner 
form of the local density of states provided that the en- 
ergy dependence of r a p is weak. The imaginary part rf^j 
is the width of the Lorentzian. In the following, we replace 
H by H , and we first evaluate ([[J]) for the case of vanish- 
ing disorder. For this we only need the projected Green 
function (^) due to the appearance of U in (|l4|). In view 
of Eq. (El) , we first determine 



(go)xy = < xx\ (E + tO-Ho) \yy > (15) 



dk e lk{x - y) g {k) 



with 



9o(k) 



dq ■ 



1 

27r J_ v E + iO + 2[cos(g) + cos(fc - q)} 
1 



(16) 



y/Acos 2 (k/2) -E 2 
The Green function at W = and U ^ is then given by 

From this and Eq. (fl4|), we obtain 

r al3 = -2U 2 ^ <Pa( x ) 9xy <Pa{y) fpiv) ■ (18) 

Inserting the plane wave eigenstates for ip a or ipp, we ex- 
actly recover the result of Ref. uM for the Breit-Wigner 
width. This shows that the diagrammatic approach of |jl4j 
is equivalent to our above approximations (replacing H by 
H and continuum limit for k). 

The generalization to the disordered case essentially 
gives rise to two modifications. Using diagrammatic per- 
turbation theory in the disorder, one can first evaluate 
the average of the Green functions ( |l5| ) and ( |l6| ) which 
amounts to the replacement E + iO — ► E + ij where 7 is 
determined by a Dyson equation. For weak disorder one 
finds 7 ~ W 2 (with eventual logarithmic corrections at 



the band center). In the following discussion, we neglect 
the effect of this small 7 which essentially regularizes ( |l6| ) 
at the singularity. The second, more important, modifi- 
cation due to finite disorder concerns the eigenfunctions 
(fa (x) in Eq. (|l8|) . Here the energy dependence of the one 
particle localization length plays an important role. To see 
this, we use the toy ansatz 



<p a (x) 



3— \&— Xot\/Lx l0l -\-i k Q 



y/Li 



(19) 



with E a = — 2cosfc Q and an energy dependent localiza- 
tion length Li jQ , L\ sin 2 (fc Q ) fl. This ansatz essentially 
corresponds to a particle moving ballistically with a well 
defined momentum inside the localization domain around 
x a . This is indeed reasonable because in one dimension 
the mean free path is of the same order as the localization 
length. However, the momentum has to be larger than its 
typical uncertainty, i. e. |fe a | > Ak ~ L^ a . Therefore the 

ansatz ( |l9| ) is valid for momenta with |sinfc Q | > 
corresponding to energies not being close to the band cen- 
ter. Inserting ( |l9| ) in ( |l8| ) and choosing x a f=a Xp, we obtain 



r, 



a/3 



-2U 1 



1 



g (k a + hp) 



Li, a + L 2 , a 1-U g (k a + kp) 
giving rise to the Breit-Wigner width (for E = 0) 
_ t/ 2 /(4Li) \cos[(k a + kp)/2}\ 



(20) 



r 



a0 sin 2 k a + sin 2 kp cos 2 [(fc Q + kp)/2] + (C//4) 2 ' 

These expressions differ by the fc-dependent localization 
length from the result of Ref. ]l4|] . Since depends 
strongly on the momenta k a and kp, we determine the 
average with respect to these momenta fl4| : 



r= 4^' dkr 



dkp5{E-E a ~Ep)/v 2 {E) . (22) 



Using Eq. (|10|), relating r with L 2 , we can estimate for 
E — > the two particle localization length as 



- 2 



L 2 « c(U) L 

2 f n 1 
c(U)=A^r dk- — 5— 
V ' ~- ./,, (sin 2 fc 1 ; " ! 

2A\U 



(sin 2 k + k 2 ) sin 2 k + (U/4) 2 

u 
1 



(23) 



irk c 



\U\s{^)^(kj + Ak cS (k c )J7(^) 



where s(x) = 1 + x 2 , A is a numerical prefactor and k c is 
a cutoff value to re gula rize the integral for small k where 
the ballistic ansatz (|l9|) is invalid. For the values k c = 0.25 
and A = 0.028 the JJ-dependence of (|23| ) fits very well the 
numerical data for the slope c(U) which can be seen in 
Fig. |5|. We can considerably simplify the somewhat lengthy 
expression ( |23| ) by neglecting the quadratic corrections U 2 
and k 2 and slightly modifying A, 



c{U) w 0.074 



1^1 



\u\ + i 



(24) 
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This approximation is numerically very accurate with an 
relative error smaller than 1 % for < \U\ < 2.0. The 
linear behavior of c(U) for small \U\ is due to a com- 
bination of the logarithmic singularity in the density of 
states at E — > and of large values of r^) ~ U if 

COs[(fc Q + kp)/2] PS ±U. 

In view of the agreement between the numerical data 
and the theoretical estimate for c(U), we conclude that 
the idea of diffusively moving particle pairs for short time 
scales finally becoming localized due to quantum interfer- 
ence P,^,^n can indeed quantitatively explain the der- 
ealization effect. For this it is important to evaluate care- 
fully the Breit-Wigncr width by taking into account the 
energy dependence of the one particle localization length. 

Despite this agreement we want to emphasize that the 
result (P3[) is nevertheless based on several qualitative ar- 
guments. Actually, the application of the relation jTo| ) 
is somewhat problematic because both the Breit-Wigner 
width and the one-particle localization length do not have 
unique values due to their energydependence. It is a priori 



not clear if the simple average (22) is really sufficient and 
accurate. Furthermore, Eqs. ([2 3D, (53) depend on the ar- 
tificial cutoff parameter k c . Theoretically, we expect that 

— 1/3 

k r ~ L 1 because of the invalidity of the ballistic ansatz 
(|19f) for small momenta |fc| < k c . Numerically, the case 
k c = 0.25 indeed corresponds to L\ ps 50 — 100 the largest 
considered L\ values. However, the resulting dependence 

7 /3 

on L\, i. e. L2 ~ clearly disagrees with the numerical 
data. We attribute this to the fact that for small momenta 
according to L\^ a ps L\ sin 2 k a the effective size of the ran- 
dom hops is strongly reduced. This feature is not properly 
taken into account in the relation . Therefore it would 
be interesting to carefully generalize this relation to the 
case where f and L\ have complicated distributions in- 
stead of unique values. 



4 Conclusion 

In this work, we have presented and applied two new accu- 
rate and efficient variants of the Greens function method 
originally introduced by von Oppen et al. |l3| to study 
the TIP localization problem ||. Our results for the TIP 
localization length L2 can be well fitted (Fig. |D by the 
functional form L2 = £i[0.5 + c(U) Li}. The behavior of 
the slope c(U) ps 0.074 \U\/(1 + \U\) is determined by 
the /7-dependence of the Breit-Wigner width r of non- 
interacting pair states. For this we presented an accurate 
estimate of r extending former work of Jacquod et at. 

We think that our results provide important additional 
evidence for the derealization effect as such. In particular, 
we find for U = 2.0 and W = 1.0 an enhancement factor 
1/2/(2 ii) pa 11.5. Our data is in qualitative agreement 
with former results of Song et al. p2fl . who directly used 
the less efficient recursive Green function technique for 
smaller system sizes (N < 200), and with very recent work 
po],p6| based on the decimation method (N < 251 and 
N < 300). In view of this the original claim of Romer et 



al. [£3| that there is no derealization effect for infinite 
system size can no longer be maintained. To understand 
the transfer matrix data on which this claim was based, 
we remind that the considered disorder value W = 3.0 
was relatively large such that L2 and L\ are nearly equal. 
Using, the Green function method one can still measure 
the enhancement because L2 > L\/2. However, in the 
transfer matrix approach there is a direct competition of 
L2 with Li |2q] . Furthermore, even for smaller disorder, 
the finite size behavior of the transmission eigenvalues is 
very subtle and one has to be careful about the finite size 
extrapolation here p2j . 

While the derealization effect is now well established, 
the situation is less clear concerning the functional de- 
pendence of L2 on L\. The formerly observed powerlaw 
L 2 ~ Lf with a ps 1.45 - 1.65 0,||,|2| was obtained 
by a fit ansatz without constant term. According to our 
above discussion (see Fig. ||) it is numerically not obvious 
to distinguish this powerlaw from the functional form we 
proposed above. Waintal et al. |^o) gave an argument in 
favor of the former with a = 1.5. This argument is based 
on the multifractal properties of the interaction induced 
coupling matrix elements in combination with an estimate 
of r using Fermi's golden rule. We believe that this anal- 
ysis is indeed important and very relevant to the prob- 
lem. Actually, our result ( pl| ) for the Breit-Wigner width 
contains a strong dependence on the initial one-particle 
states due to partial momentum conservation. This leads 
to strong fluctuations of r which are presumably directly 
related to the multifractal statistics of the coupling matrix 
elements. However, the analytical calculations of Ref. [ fu] 
and of section ^| clearly show that the simple application 
of Fermi's golden rule is not sufficient and the numerical 
data do not show the corresponding behavior L2 ~ U 2 . 
To understand these issues in more detail further work is 
necessary. 

We emp hasize that our numerical data and the es- 
timate (G3p are valid for moderate interaction strengths 
\U\ < 2.0. For larger values of U one expects that L2 will 
be reduced due to the particular projector structure of 
the interaction operator U. Waintal et al. ]2(l presented 
a duality transformation mapping the case of \U\ 1 to 
a similar problem with |{7| <C 1 and a different reference 
basis. According to this Li2(U) should obey the duality 
relation L 2 (U) ps L 2 {V2A/U) @. 

Finally, we mention that the numerical trick to eval- 
uate efficiently the matrix go via Eq. (Q) also works in 
higher dimensions, even though the gain is less spectacu- 
lar. In d dimensions and a system of total size (volume) 
N one can calculate G {1 \E) by the recursive Green func- 
tion method which provides an algorithm to evaluate (0) 
with 7V 3+ ( d ~ 1 )/ d operations. In particular the case of two 
dimensions is important due to recent claims of Ortuno 
et al. and Cuevas ]33| for a derealization transition for 
two interacting particles in d — 2. We think it is necessary 
to consider larger system sizes as in Ref. |33| in order to 
decide whether there is a real transition or a very strong 
derealization with a finite but very large two-particle lo- 
calization length as it was argued by Shepelyansky M. 
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